GW method with the self-consistent Sternheimer equation 
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We propose a novel approach to quasiparticle GW calculations which does not require the com- 
putation of unoccupied electronic states. In our approach the screened Coulomb interaction is 
evaluated by solving self-consistent linear-response Sternheimer equations, and the noninteracting 
Green's function is evaluated by solving inhomogeneous linear systems. The frequency-dependence 
of the screened Coulomb interaction is explicitly taken into account. In order to avoid the singular- 
ities of the screened Coulomb interaction the calculations are performed along the imaginary axis, 
and the results are analytically continued to the real axis through Pade approximants. As a proof of 
concept we implemented the proposed methodology within the empirical pseudopotential formalism 
and we validated our implementation using silicon as a test case. We examine the advantages and 
limitations of our method and describe promising future directions. 

PACS numbers: 71.15.-m, 71.15. Qe 



I. INTRODUCTION 

During the past two and a half decades the GW 
method 1,2 for the study of electron quasiparticle exci- 
tations has had a number of successes and witnessed 
significant growth of interest within the computational 
electronic structure community. The GW method is cur- 
rently being used for predicting electron quasi-particle 
excitation spectra as well as optical spectra in a vari- 
ety of materials ranging from bulk solids to nanostruc- 
tures and organic systems. The GW method is also of 
widespread use as a starting point for Bethe-Salpeter cal- 
culations of two-particle neutral excitations. 3-8 Current 
implementations find many diverse applications, includ- 
ing among others the calculation of the optical response 
of nanostructures, 9 quantum transport in nanoscale 
junctions, 10 pump-probe spectroscopy, 11 angle-resolved 
photoemission spectroscopy, 12 and strongly-correlated 
systems. 13 

Current trends in the development of improved com- 
putational approaches for quasiparticle excitations based 
on the GW method include the refinement of the initial 
guess for the non-interacting Green's function and for the 
polarization operator, 14,15 the inclusion of approximate 
vertex corrections or higher-order self-energy diagrams, 16 
and the description of the frequency-dependent dielectric 
response beyond the original generalized plasmon-pole 
approximation. 17 ' 18 Detailed reviews of past and current 
developments in GW techniques can be found in Refs. 
8,14,19-21. 

The majority of current GW implementations ob- 
tain the screened Coulomb interaction W and the non- 
interacting Green's function G using a perturbative ex- 
pansion over the Kohn-Sham cigenstates (cf. Sec. II A be- 
low). Such expansion requires the calculation of both oc- 
cupied and unoccupied electronic states, as well as their 
associated optical matrix elements. 2 A common bottle- 



neck of this approach is that the convergence of the quasi- 
particle excitation energies with the number of unoccu- 
pied states is rather slow. 22 This difficulty is particu- 
larly relevant when calculating the absolute values of the 
quasiparticle excitation energies. 23 Several avenues have 
been explored so far in order to circumvent this bottle- 
neck and to perform GW calculations by employing only 
occupied electronic states, 24-27 or a small number of un- 
occupied states. 23 

The main aim of the present work is to demonstrate 
the feasibility of GW calculations entirely based on oc- 
cupied states only. 28 In practice we adopt the princi- 
ples of density-functional perturbation theory (DFPT) 
to determine (i) the frequency-dependent screened 
Coulomb interaction by directly solving self-consistent 
linear response Sternheimer equations, and (ii) the non- 
interacting Green's function by solving inhomogeneous 
linear systems. The main advantage of the proposed 
method is that it does not require the computation of un- 
occupied electronic states. In addition, we demonstrate 
the possibility of fast evaluations of the frequency de- 
pendence of the screened Coulomb interaction based on 
multishift linear-system solvers. 29 As a proof of concept 
we have implemented our method within a planewaves 
empirical pseudopotential scheme, 30 and validated it by 
comparing with previous work for the prototypical test 
case of silicon. 

The use of the Sternheimer equation for calculating the 
polarizability in the random-phase approximation (RPA) 
or the inverse dielectric matrix has already been discussed 
in Refs. 31 and 32, respectively, within the framework 
of a non-perturbative supercell approach. After the in- 
troduction of DFPT in the context of lattice-dynamical 
calculations, 33 the authors of Ref. 24 proposed the use 
of the non self-consistent Sternheimer method for the 
calculation of the dielectric matrix. The elimination 
of unoccupied electronic states in the evaluation of the 
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screened Coulomb interaction has also been proposed re- 
cently within the framework of a Wannicr-like represen- 
tation of the polarization propagator and the Lanczos 
recursion method. 25,26 

This manuscript is organized as follows. In Sec. II we 
describe how the self-consistent Sternheimer formalism 
can be adapted to perform GW calculations. In partic- 
ular, we outline the procedure to obtain the screened 
Coulomb interaction in Sec. II A, the non-interacting 
Green's function in Sec. II B, and the GW self-energy 
in Sec. II C. In Sec. Ill we specialize to a planewave ba- 
sis set representation and derive the key equations for 
the case of Bloch electrons. Sections III A, IIIB, IIIC 
parallel the corresponding sections in the general theory 
part, respectively. In Sec. IV we critically analyze the ad- 
vantages and limitations of the present approach with an 
emphasis on the scaling of the calculations with system 
size. In Sec. V we report the results of our test calcula- 
tions for silicon and compare with previous calculations 
in the literature. Specifically, we presents results for the 
direct and inverse dielectric matrix (Sec. VA), for the 
analytic continuation of the dielectric matrix using Pade 
approximants (Sec. VB), for the self-energy (Sec. VC), 
and for the spectral function (Sec. VD). In Sec. VI we 
discuss possible future developments of our method and 
discuss our conclusions. The Appendices provide tech- 
nical details on some numerical algorithms adopted in 
this work, in particular the preconditioned complex bi- 
conjugate gradient method (Appendix A), the analysis 
of the conditioning of the Sternheimer equations (Ap- 
pendix B), the analytic continuation using Pade approx- 
imants (Appendix C), and the use of multishift methods 
for the simultaneous calculation of the polarization at 
multiple frequencies (Appendix D). 



II. GENERAL THEORY 

A. Screened Coulomb interaction 

In this section we describe how to exploit the Stern- 
heimer scheme within density-functional perturbation 
theory in order to calculate the screened Coulomb in- 
teraction W(r, r';ui) (where r, r' arc the space vari- 
ables and 10 is the excitation frequency). While the 
use of the Sternheimer approach in DFPT was origi- 
nally developed bearing in mind the Kohn-Sham effective 
Hamiltonian, we note that the present procedure applies 
without restrictions also to post-DFT methods such as 
the LDA+U method, 34 hybrid functionals, 35 and exact 
exchange. 36 We assume Rydberg atomic units through- 
out this manuscript. The Hedin's equation which defines 
the screened Coulomb interaction reads: 20 

W(r,r';u;) = v(r,r') + J dr" W(r, r"; w) 

x f dv"'P{v",v"'^)v{v"',v'), (1) 



where v(r, r') = e 2 /|r — r'| denotes the bare Coulomb 
interaction and P(r, r';u>) the irreducible polarization 
propagator. As Eq. (1) is a self-consistent Dyson equa- 
tion for the screened Coulomb interaction, it should be 
possible to solve it recursively in the spirit of density- 
functional perturbation theory. For simplicity, we here 
specialize to the case of the random-phase approxima- 
tion (RPA) for the polarization propagator. The gener- 
alization of this procedure to include exchange and cor- 
relation effects can be performed without difficulties (cf. 
Sec. II A 1). Within the random-phase approximation the 
polarization propagator can be written as: 20 

P(r,r';c) = 2 ]T /" = fm 4 (^(^(0^,(0, 

nm 

(2) 

where tp n (r) indicates an electronic eigenstate of the 
single-particle Hamiltonian with energy eigenvalue e n 
and occupation number /„. In the following we assume 
that the tl> n (r) are Kohn-Sham eigenstates for definitc- 
ness. In Eq. (2) the summation indices m and n run over 
both occupied and unoccupied electronic states, and the 
factor of 2 accounts for the spin degeneracy. 20 Although 
the expression for the RPA polarization Eq. (2) has been 
derived for real frequencies in Ref. 20, it is possible to 
continue the polarization throughout the complex plane 
by using Eq. (2) as a definition outside of the real axis. 

Our goal is to rewrite Eqs. (1) and (2) by avoid- 
ing explicit summations over the unoccupied electronic 
states. For this purpose it is convenient to regard the 
screened Coulomb interaction VF(r,r';u>) as a function 
of the second space variable r', whilst the first space 
variable r and the frequency u> are kept as parameters: 
AVj FjW ](0 = W(r,r';u>). If the system under consid- 
eration is subject to the perturbation AVj r)W ](r'), then 
within the RPA the first-order variation of the single- 
particle density matrix Anr riW i reads 

An M = 2j2r v ^r v[ rM- ( 3 ) 

va 

In Eq. (3) the index v stands for "valence" and runs over 
the occupied electronic states only, the factor of 2 is for 
the spin degeneracy, and the superscript a = ± refer to 
the positive and negative frequency components of the 
induced charge. The first-order variations of the occupied 
wavefunctions A^£Q r ^ can be determined by solving the 
following two Sternheimer equations: 

(H — s v ± u)M% M = -(1 - Pocc)AV [ria)] ^, (4) 

where H is the effective single-particle Hamiltonian and 
-Pocc = J2 V \i } v){i ) v\ is the projector on the occupied 
states manifold. In the particular case of vanishing fre- 
quency (uj = 0) the a = ± variations of the wavefunctions 
do coincide, and the standard DFPT equations are recov- 
ered. The screening Hartree potential associated with the 
induced charge An[ FjW ] is calculated as usual through 

AV£ w] (0 = / <ir"An M (r>(r",r'), (5) 
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and finally the screened Coulomb interaction in the RPA 
is obtained as 

W(r, r'; u) = AV [rM (r') = v(v, r') + AV£ u] (r'). (6) 

It is tedious but otherwise straightforward to verify 
that Eqs. (3)-(6) are equivalent to the original Eqs. (1)- 
(2). The only assumptions made in our derivation are 
that time-reversal symmetry applies, and that the system 
under consideration has a finite energy gap for electronic 
excitations. The assumption of time-reversal symmetry 
is not essential and is mainly used to obtain a compact ex- 
pression for the a = ± wavefunction perturbations. The 
assumption of finite energy gap can be relaxed by using 
the extension of DFPT to metallic systems developed in 
Ref. 37. 

There is a simple and intuitive physical meaning as- 
sociated with the calculation scheme outlined above. To 
see this we consider an external test charge introduced 
in the system at the point r. This charge generates 
a bare Coulomb potential v(r,r'), and the system re- 
sponds to such perturbation by generating the induced 
charge An[ rjW ](r') and the associated screening poten- 
tial AV^ w j(r'). The sum of the external perturbation 
v(r, r') and the screening potential AV^ ( r ') yields the 
screened Coulomb interaction W(r,r';u>) at the point r' 
within the RPA. 

The linear systems in Eq. (4) must be solved self- 
consistently. For this purpose we begin by initializing 
the screened Coulomb interaction W using the bare in- 
teraction v. We then calculate the linear variations of 
the wavefunctions Atp°. Using the calculated linear vari- 
ations we update the induced charge density An and the 
associated screening potential AF H . This allow us to 
generate an improved estimate of the screened Coulomb 
interaction W. We cycle through these steps until con- 
vergence of the screened Coulomb interaction is achieved. 
Equations (3), (4) can be regarded as the generalization of 
the self-consistent Sternheimer equations used for lattice- 
dynamical calculations 28 to finite-frequency test-charge 
perturbations. 

In practical calculations we solve Eq. (4) along the 
imaginary frequency axis in order to avoid the null eigen- 
values of the operator H — e v ± u>, and then we perform 
the analytic continuation of the screened Coulomb inter- 
action to real frequencies (cf. Appendixes A-C). In the 
special case of uj = it is convenient to modify the linear 
operator on the left-hand side of Eq. (4) by adding the 
projector on the occupied states manifold P c C : 

(H-e v + aP occ )A^ v[rfi] = -(1 - P occ )AVj r;0] W, (7) 

with a set to twice the occupied bandwith. This ex- 
tra term does not affect the solutions A-i/^^o] which are 
linear combinations of unoccupied electronic states. At 
the same time, the extra term has the effect of shifting 
away from zero the null eigenvalues of the linear operator 
H — e v thereby making it non singular. This strategy is 
common practice in DFPT implementations, 28 ' 38 and is 
discussed in greater detail in Appendix B. 



1. Vertex correction 

Within the scheme outlined here it is rather straight- 
forward to introduce an approximate vertex correction 
to the GW self-energy along the lines of Refs. 2,39. This 
correction results from setting the self-energy in the first 
iteration of Hcdin's equations to the DFT exchange- 
correlation (XC) potential, So(r, r';cu) — S(r, r')V xc (r). 
Within the present scheme this correction is simply 
obtained by including the variation of the exchange- 
correlation potential in the self-consistent potential used 
in Eq. (4): 



(fl"-e„±w)A^ [P)U] =-(l-P cc) AV [rM +K xc An [r ^ 



(8) 

K xc = SV xc /Sn being the functional derivative of the 
XC potential with respect to the density. The screened 
Coulomb interaction is still to be calculated through 
Eq. (6). This approach has been called "GW + K xc " 
approximation in Ref. 39 due to the inclusion of the XC 
contribution in the screening of the test charge. That 
the inclusion of the XC term in the self-consistent in- 
duced potential leads to the GW + K xc approximation 
can easily be seen as follows. We combine Eqs. (3), (8) to 
yield the induced charge density (we use symbolic oper- 
ator notation for clarity): 



An = v[l-P(v + K xc )]- 1 P. 



(9) 



Then, we substitute this result in the definition of the 
screened Coulomb interaction Eqs. (5), (6) to find 

W = v{l + v[l-P(v + K xc )]- 1 P}. (10) 

The last equation yields precisely the screened Coulomb 
interaction in the GW + K xc approximation. 2 ' 39 The dif- 
ference between this approach and the standard GW ap- 
proximation is that in this case the screening charge is 
calculated for an electron, while in the GW-RPA approx- 
imation the screening is calculated for a test charge. It 
is worth pointing out that in standard implementations 
of DFPT the XC term is already included in the varia- 
tion of the self-consistent potential, 28 therefore the use 
of the GW + K xc approximation would not require any 
additional computational developments if the present ap- 
proach was to be implemented on top of existing DFT 
software. 



2. Non self- consistent calculation of the dielectric matrix 

An alternative approach to the calculation of the 
screened Coulomb interaction using the self-consistent 
Sternheimer equation consists in solving Eq. (4) non 
sclf-consistently. For this purpose we can replace the 
self-consistent perturbation AVj r , w ](r') in the right-hand 
side of Eq. (4) by the bare Coulomb potential f[ r ](r') = 
v(r, r') as follows: 

(H-S V ± W)A^;J = -(1 - PoccMr]^, (11) 
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and we can solve this Stcrnheimcr equation with the 
known term on the right-hand side kept fixed. By con- 
structing the non self-consistent induced charge density 
An^ s w , as in Eq. (3) we then obtain the dielectric matrix 
e(r,r';w): 

e(r,r';u,) = S(r,r>)-Anf r % ] (r'). (12) 

It is straightforward to check that this procedure cor- 
rectly leads to the RPA dielectric matrix. 40 The differ- 
ence between this approach and the self-consistent calcu- 
lation described in Sec. II A is that here we also need to 
invert the dielectric matrix obtained through Eq. (12) in 
order to calculate the screened Coulomb interaction. 

This non self-consistent procedure was first proposed in 
Ref. 24. One additional step that we make in the present 
work is to notice that Eq. (11) constitutes a shifted linear 
system, i.e. a system where the linear operator H — s v ±lo 
differs from the "seed" operator H—e v only by a constant 
shift ±loI (I being the identity operator). In this case we 
can take advantage of the multishift linear system solver 
of Ref. 29 to determine A^^'^j for every frequency lo 
at the computational cost of one single calculation for 
the seed system (H - e^A^^o] = _ (1 ~ Pocc)v[ r ]ipv 
This procedure is extremely advantageous as it makes it 
possible to calculate the entire frequency-dependence by 
performing one single iterative minimization. The tech- 
nical implementation of this procedure is described in 
Appendix D. 

B. Green's function 



In the above derivation we assumed again time-reversal 
symmetry, we used the Lorentzian representation of the 
Dirac's delta function for small rj [ttS(x) = r]/(x 2 + rj 2 )], 
and we defined e~ = e n — ir]. The component G A of 
the Green's function is obviously analytic in the upper 
half of the complex energy plane as its poles lie below 
the real axis. The non-analytic component G N vanishes 
whenever the frequency lo is above the chemical poten- 
tial. For frequencies lo below the chemical potential, the 
non-analytic component introduces the poles associated 
with the occupied electronic states. The partitioning of 
Eqs. (16), (17) closely reflects the analytic structure of 
the non-interacting Green's function. A detailed discus- 
sion of this aspect can be found in Ref. 20. The two 
components of the Green's function in Eq. (16) are asso- 
ciated with the Coulomb hole (COH) and the screened 
exchange (SEX) terms of the self-energy, £ COH = G A W 
and E SEX = G N VK, respectively. 2 

The computation of the non-analytic component G N 
of the Green's function in Eq. (17) is straightforward 
once the occupied electronic eigenstates have been de- 
termined. In order to calculate the analytic component 
G A it is convenient to proceed as in the case of the 
screened Coulomb interaction, by regarding G A (r, v';lo) 
as a parametric function of the the first space variable 
and of the frequency: G^^j(r') = G A (r, r'; lo). If we 

apply the operator H — lo + to both sides of Eq. (16), 
with lo + = lo + ir], and we use the completeness rela- 
tion £[ r ](r') = 6(r,r') = E n ^n( r )C( r ')i tncn wc nnd 
immediately: 



The calculation of the Green's function can efficiently 
be performed by adopting a strategy similar to the Stern- 
heimer approach described in Sec. II A. We introduce the 
noninteracting Green's function following Ref. 20: 



G(r,r» = £ 



Vqr)C(r') 
lo- e n - ir]„ 



(13) 



where the sum extends over occupied as well as unoccu- 
pied electronic states. The real infinitesimal rj n is positive 
(Vn = ff) for occupied states and negative {j] n = —rj) for 
unoccupied states. 2,20 ' 41 We now split the sum in Eq. (13) 
into occupied (v) and unoccupied (c) states: 

G (r,r>).V«) + V™, (14) 
^-^ lo — e v — in ^ lo — e c + ir> 

V ' c ' 

and we add and subtract ^2 V rp v rp*/(Lo—s v +ir]) to obtain: 
G(r, r'; lo) - G A (r, r'; lo) + G N (r, r'; lo), (15) 

with 

^(r)V-n(r') 



G a (v,v';lo) = 



LO - E r , 



(16) 



G N (r,r';^) = 2™ £ 5(lo - e^r)^)- (17) 



{H-lo+)G{ 



(18) 



As expected, we can determine the analytic part of the 
Green's function by directly solving a linear system. As 
Eq. (18) does not explicitly require unoccupied elec- 
tronic states, this procedure mimics the Sternheimer ap- 
proach for the screened Coulomb interaction outlined in 
Sec. II A, albeit without the self-consistency requirement. 

The procedure described here is especially advanta- 
geous because Eq. (18) constitutes a shifted linear sys- 
tem, in the same way as Eq. (11). Also in this case we 
exploit the multishift method of Ref. 29 to determine 
Gj^; w j for every frequency lo at the computational cost of 

one single calculation for the seed system HG^ j = — <5[ r ] 
(cf. Appendix D). 

The presence of the infinitesimal it] in lo + = to + ir] 
guarantees that the linear operator H — lo + in Eq. (41) 
is never singular. This operator can nonetheless become 
ill-conditioned, hence the use of appropriate precondi- 
tioners may become necessary. We discuss this aspect in 
Appendix B. 
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C. Self-energy 



The electron self-energy in the GW approximation is: 2 



S(r,r'; t 



du'G(r, t',oj + u')W(r, r', u')e 



-iSuj' 



(19) 

where 8 is a positive infinitesimal. At large frequen- 
cies the Green's function decays as uj^ 1 and the screened 
Coulomb interaction tends to the frequency-independent 
bare Coulomb interaction v. As a consequence, the in- 
tegrand in Eq. (19) decays as o> _1 and the integration 
requires some care. 

It is convenient to split the self-energy into an ex- 
change contribution S cx = Gv and a Coulomb term 
S c = G(W — v). 18 It is easy to check that the integrand 
in the Coulmb term decays as u>~ 2 at large frequencies, 
therefore the integral is well behaved and the integration 
can be performed by using a numerical cutoff \u>'\ < u>c 
in Eq. (19). A detailed analysis of the analytic properties 
of the Coulomb term S c (w) shows that it must decay as 
lu^ 1 at large frequencis, and that the use of the cutoff 
luq in the integration introduces an error of the order 
of ojp/uje, where oj p denotes the characteristic plasmon 
frequency of the system. 

In order to integrate the exchange term we observe 
that S ox = G A v + G N v and that the poles of G A lie 
in the lower half of the complex plane, hence the inte- 
gration of the term G A v yields a vanishing contribution. 
On the other hand, the integration of G N v yields a con- 
stant (frequency-independent) term. 42 In summary, we 
perform the frequency integration in Eq. (19) by evaluat- 
ing numerically the Coulomb term using an energy cutoff, 
and by integrating analitycally the exchange term: 



E(r,r';a;) = S c (r,r';u;) + S ex (r,r'), 



(20) 



with 



S c (r,r';w)=— / dw'G?(r, r'; u+u') W{r, r', J) - v(r, r' 

2?r "'— c 

(21) 

and 

E cx (r,r') = -^^(r)^(r')»(r,r'). (22) 



III. IMPLEMENTATION IN A BASIS OF 
PLANEWAVES 

We here describe our planewaves implementation of 
the scheme developed in Sec. II. The choice of a 
planewaves representation was motivated by the need 
for making contact with existing literature on dielectric 
matrices, 40 ' 43 " 45 and by the availability of DFPT soft- 
ware for lattice-dynamical calculations 38 which was used 
as a reference for our implementation. 



We adopt the following conventions for the transfor- 
mation from real to reciprocal space. The wavefunctions 
transform as usual according to: 

^ k (r') = JL e*( k+G ')' r ' M „ k (G'), (23) 

G' 

with fi the volume of the unit cell and k the Bloch 
wavevector. The bare Coulomb interaction transforms 
according to 

v(r,r') = ^$>(q+G)c^+ G H r '- r ), (24) 

q qG 

where q is also a Bloch wavevector and _/V q is the number 
of such wavevectors in our discrctized Brillouin zone. In 
Eq. (24) we have v(q + G) = 47re 2 /|q + G| 2 . The latter 
expression for v(q + G) is arrived at by replacing the 
integration J dr exp(iq • r) / |r | over the crystal volume by 
an integration over all space. This choice corresponds to 
assuming that we can rely on a very fine sampling of the 
Brillouin zone. Had we performed the integration on a 
sphere with a radius R c defined by the crystal volume 
(4/37ri? 3 = NqQ), then we would have obtained 



t*(q + G) 



47re 2 



|q + G| 2 



(l-cos|q + G|i? c ), (25) 



which corresponds to the truncated Coulomb potential 
introduced in Rcf. 46. We will come back to this aspect 
in Sec. Ill A. The screened Coulomb interaction and the 
non interacting Green's function transform according to 



W{v,v'-uj) = 
and 

G(r,r';u) 



1 



1 



e -*(q+G)-r WGG , (q;a;)e i(q+G')-r' 5 

(26) 



kGG' 



e -i(k + G).r GGG/(k;w y(k + G')-r' ; 

(27) 

with similar expressions for G A , G N . We note that the 
sign convention adopted here in the Fourier transforms 
[e.g. exp(+?q • r') in the rhs of Eq. (26)] is necessary to 
obtain the compact expression Eq. (35) below for the in- 
duced charge, and is opposite to the convention adopted 
in Ref. 2. Before proceeding it is also convenient to intro- 
duce the "right-sided" inverse dielectric matrix through 



(28) 



By adopting the same convention for the inverse dielectric 
matrix as for the screened Coulomb interaction the above 
equation can be rewritten as: 



FT GG '(q;w) = v(q + G)e GG ,(q; u). 



(29) 



We note that our Eq. (29) is slightly different from the 
standard expression [e.g. Eq. (22) of Ref. 2], due to our 
choice of using the right-sided inverse dielectric matrix. 
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A. Screened Coulomb interaction 

In order to rewrite Eqs. (3)-(6) in the Bloch represen- 
tation and in reciprocal space we proceed as follows. We 
first write the linear systems Eq. (4) by relabeling the 
wavefunctions tjj v as Bloch states tp^: 



k[i\w] 



-(l-P occ )AV[ r;W] ^ k . (30) 



The linear variations of the wavefunctions can be ex- 
panded in terms of Bloch waves as follows: 

A,/, ff , , — 1 Av a , ,p i ( k +q)- r ' P ~ i ( c '+ G )' r 

JVc l" qG 

(31) 

where Aw^ k j q G , is cell-periodic in r'. From the linear 
variations of the wavefunctions we construct the induced 
charge density using Eq. (3): 



An 



■E^Ac k 



1 ■-' ^ YviK -^ Yv ^[ r M' 

i?ktx 



(32) 



Here the factor iVk takes into account the normalization 
of the Bloch states in the unit cell [the wavefunctions ip v 
in Eq. (4) are normalized in the whole crystal]. Next 
we expand the screened Coulomb interaction in terms of 
Bloch waves: 

A^(r') = ^^A^ G , w] (r')e-^ G )-V q - r ', 

JVq " qG 

(33) 

where Ai>[ qjG;W ] is cell-periodic in r'. If we now place Eqs. 
(33) and (31) into Eq. (30) we discover that the compo- 
nent Aw[ q G w ] of the perturbing potential corresponding 
to the Bloch wave exp(— iq • r) couples only to the vari- 
ations of the wavefunctions corresponding to the Bloch 
wave exp[z(k + q) • r']. As a result the linear system Eq. 
(4) becomes: 

(i?k+ q -et,k±w)Aii± k[qG w] = -(l-P k + q )Av [q!G ^ ]Mt , k , 

(34) 

where H k = e -* kr 7ie 4kr and P k cc = \u vk )(u vk \. The 
induced charge density associated with the Bloch wave 
exp(— iq • r) now reads 

2 x - 

An [q , G ,a,] = Z^ u "kA< k[qiG>w] . (35) 



This result is very similar to the case of standard 
DFPT. 28 The main difference is that in the present case 
the translational invariance of the screened Coulomb in- 
teraction induces a coupling between the perturbation 
with Bloch wave exp(— «q • r) in the variable r and its 
induced response with Bloch wave exp( +iq • r') in the 
variable r'. To conclude our derivation, we rewrite the 
screened Coulomb interaction Eq. (6) after expanding the 
cell-periodic function An[ q>GiW ] (r') in planewaves: 

WW(q;w) - [<5 GG < + An [q , GjW] (G')Mq + G')- (36) 



In practical calculations we proceed as follows: we 
first initialize the perturbation in the linear systems us- 
ing A^ Gj ,(r') = v(q + G)exp(iG • r'). The solu- 
tion of the linear systems yields the change in the wave- 
functions, which are then used to construct the induced 
charge, the induced Hartree potential, and the updated 
screened Coulomb interaction. We repeat this procedure 
by starting with the updated screened Coulomb interac- 
tion until convergence is achieved. At convergence the 
self-consistent perturbing potential yields the screened 
Coulomb interaction W GG '(q;w). The calculation must 
be repeated for every perturbation, i.e. for each set of 
parameters [q, G,ui}. At the end of this procedure it 
is straightforward to obtain the inverse dielectric matrix 
e GG ,(q;w) using Eqs. (36) and (29). Alternatively, it 
is also possible to scale the initial perturbation and use 
Ad^q w j = exp(zG • r') to obtain the inverse dielectric 
matrix at the end of the self-consistent procedure [indeed 
Eq. (34) is a linear system]. 

The scheme developed here allows us to calculate one 
row (in G') of the inverse dielectric matrix e GG ,(q;o;) 
by determining the linear response to the perturbation 
exp[— i(q+ G) • r]. This idea has been discussed already 
in Ref. 32 in the framework of nonperturbative methods 
based on supercell calculations. 



1. Singularities in the inverse dielectric matrix and the 
screened Coulomb interaction 

In order to avoid the singular behavior of the wings of 
the inverse dielectric matrix in the long wavelength limit 
(|q| — >0) it is convenient to work with the symmetrized 
inverse dielectric matrix defined as follows: 44 



— i < \ —1/ \\ < l 

e GG'(q; w ) = e GG'(q; w )ir: 



G'l 



q+G 



(37) 



Unlike its unsymmetrized counterpart e GG , (q;w), the 
wings of e GG / (q; u>) have finite limits at long wavelengths. 
The screened Coulomb interaction of Eq. (29) is now 
rewritten in symmetrized form as 



W GG >{q:uj) = 



Aire 1 



e GG ,(q;w). (38) 



|q + G||q+G 



77 GG 



While the symmetrized inverse dielectric matrix has fi- 
nite limits at long wavelengths, the screened Coulomb in- 
teraction still presents a divergence corresponding to the 
long-range tail of the Coulomb potential in real space. 
This divergence requires special handling when perform- 
ing the Brillouin zone integration to calculate the GW 
self-energy. 2 We here overcome this difficulty following 
the prescription of Ref. 46. For this purpose we replace 
the bare Coulomb potential v(r, r') by the truncated po- 
tential v t (r, r') = u(r,r')[l-e(|r-r'|-i? c )], Q(x) being 
the Heaviside step function. The truncation radius is 
defined as in Sec. Ill A. Using this truncated Coulomb 
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potential, the final expression for the screened Coulomb 
interaction in reciprocal space becomes: 

W GG ,( q ^) = W ^U*")- (39) 

In the long- wavelength limit q^O the head of the trun- 
cated screened Coulomb interaction (G = G' = 0) tends 
to the finite limit 27re 2 i? 2 Eqq 1 (q — > 0;o>) and the singu- 
lar behavior is removed. Optimized truncation strategies 
have been developed for non-isotropic materials and sys- 
tems with reduced dimensionality. 22 

B. Green's function 

We now specialize Eqs. (13)-(18) to the case of a 
planewaves basis and the Bloch representation. We start 
by rewriting Eq. (13) after relabeling the electronic states 
ip n as Bloch states -0nk and taking into account the nor- 
malization, as already done in Sec. Ill A. Next we ex- 
pand the Green's function in terms of the Bloch waves 
exp[— i(k + G) • r] and exp(ik • r'): 

Gtr M (r') = ^£sfk,G,.](r0e- <(k+G) -V k - r \ (40) 

k kG 

with G w ](r') cell-periodic in r'. An analogous expan- 
sion holds for the non-analytic component G N . Equa- 
tions (16), (17) are now rewritten as: 

(Hu-u + )g^ GM (G') = -5 GG ,, (41) 
.9[k,G,.](G') - 27ri£%-^ k K k (G)«„ k (G'). (42) 

V 

In deriving Eqs. (41), (42) we made use once again of 
time-reversal symmetry, yielding u* k (G) = u Vt -k(— G). 
Similarly to the case of the screened Coulomb interac- 
tion, by solving the linear system in Eq. (41) for a set of 
parameters [k, G,w] we obtain an entire row G' of the 
analytic component of the Green's function g^ G w j(G'). 

C. Self-energy 

The GW self-energy in Eq. (19) is calculated in 
real space after performing the Fourier transforms of 
W^GG'(q;w) and C?GG'(k;u;). The result is then trans- 
formed back in reciprocal space to obtain Sgg'C 4 ;^). 
The evaluation of the matrix elements of the self-energy 
in the basis of Kohn-Sham eigenstates is performed in 
reciprocal space. Since the planewaves cutoff required 
to describe the inverse dielectric matrix and the self- 
energy is typically much smaller than the cutoff used in 
density-functional calculations, 2 the procedure described 
here does not require an excessive computational effort 
and accounts for only a fraction of the total computation 
time. 



IV. SCALING PROPERTIES 

In this section we analyze the computational complex- 
ity of the algorithms proposed in Sec. Ill, by focusing 
on our planewaves implementation. Without loss of gen- 
erality we consider a T point sampling of the Brillouin 
zone and we leave aside the frequency-dependence. We 
assume that the Kohn-Sham electronic wavefunctions are 
expanded in a basis of planewaves with a kinetic en- 
ergy cutoff E™* t , corresponding to N G plane waves. In 
the simplest case of norm-conserving pseudopotential ap- 
proaches the electronic charge density is described using 
a basis set with a cutoff = 4E™* t , and the corre- 

sponding numbers of basis functions and real-space grid 
points are N G n and A den , respectively. The screened 
Coulomb interaction and the Green's function are de- 
scribed by a smaller cutoff E* ut and N G planewaves. The 
self-energy is expanded in a planewaves basis with cut- 
toff E^ t = 4E£ ut , and we denote by E the number of 
real-space grid points associated with this basis. 



1. Screened Coulomb interaction 

Equation (34) needs to be solved for each one of the N G 
planewave perturbations and the N v occupied electronic 
states. For the solution of Eq. (34) we adopt the com- 
plex bi-orthogonal conjugate gradient method of Ref. 47 
(cBiCG), as described in Appendix A. Each solution 
of Eq. (34) requires two cBiCG minimizations (for ±w), 
and each cBiCG minimization consists of two conjugate 
gradients (CG) sequences. The most time-consuming op- 
eration in each CG step is the application of the Hamil- 
tonian to the previous search direction, and in particular 
the Fourier transform of the wavefunctions to real-space 
and back for evaluating the product with the local poten- 
tial. Fast-Fourier-transform (FFT) algorithms allow us 
to perform these calculation in N^ n T = 4A r don logA r den 
floating point operations. 48 If in average the CG mini- 
mization requires Nqg steps and the self-consistency loop 
requires Ascf iterations, then the total cost of the en- 
tire calculation corresponds to a number of floating point 
operations 

^!o G P f = 8N CG N SCF N G N v N& n T} (43) 

where SGW stands for "Sternheimcr GW". As N G , N v , 
and A den scale linearly with the size of the system as 
measured by the number of atoms iV at , the overall scaling 
of this procedure is A at logA at . 

For comparison it is useful to consider the scaling of 
standard GW calculations based on the expansion over 
unoccupied states (hereafter referred to as the "HL" 
method). 2 The calculation of the irreducible RPA po- 
larization requires the evaluation of the optical matrix 
elements between each of the N v occupied states and 
each of the N c unoccupied states. These matrix elements 
are typically computed by using Fourier transforms of 
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ip*(r)ip v (r), therefore this procedure requires essentially 
N V N C Fourier transforms from real- to reciprocal-space. 
Each Fourier transform is performed on the real-space 
grid for the density with N^ en grid points, therefore the 
total cost of the standard method corresponds to a num- 
ber of floating point operations 

A fl H D L ps = N c N v N« c F n T . (44) 

Even in this case therefore the overall scaling is 
A a 3 t logA at . 

Since the method of Ref. 2 calculates the dielectric 
matrix and then performs a matrix inversion, in order to 
compare the prefactors in Eqs. (43) and (44) we consider 
the non self-consistent calculation of the dielectric matrix 
as described in Sec. II A 2 [Ascf = 1 in Eq. (43)], and we 
restrict ourselves to the calculation of the static dielectric 
matrix. In this case only one calculation of Eq. (4) is 
required instead of two for ±w, and the two CG sequences 
of the cBiCG algorithm do coincide. As a result, a factor 
4 drops out of the prefactor in Eq. (43). If we assume 
for definitcness a perfectly well-conditioned linear system 
(condition number k = 1), and express the number of 
CG iterations required to achieve convergence through 
Eq. (Bl), then the ratio between the complexity of the 
SGW approach in a planewaves implementation and the 
standard approach becomes 

^ S o G pf/<o L ps = NyN c log(2/e), (45) 

where e is the relative accuracy of the results. As an ex- 
ample, for a relative accuracy of e = 10~ 5 we find this 
ratio to be ~ 12Aq/A c . In the case of silicon, using a 
typical cutoff E* ut = 10 Ry we obtain = 137, there- 
fore the SGW approach becomes convenient when more 
than ~1650 unoccupied states are used in the standard 
approach. This is rarely the case as most calculations 
reported to date use only a few hundreds of unoccupied 
electronic states. Of course the accuracy of the stan- 
dard sum-over-states expression is difficult to quantify, 
and probably a convergence on 5 significant digits is not 
warranted by a few hundreds of electronic states. 

Our estimate suggests that the planewaves implemen- 
tation of our method can be as expensive as the stan- 
dard approach. It should be noted, however, that our 
scheme has the advantage of providing the whole self- 
energy E(r, r';oj), while the standard approach typically 
provides the matrix elements of the self-energy on a small 
subset of states of the order of N v . Therefore if we were 
to perform a comparison based on the same amount of 
output information, we should use N C N^ in Eq. (44) in- 
stead of N C N V . In this case Eq. (45) would change into 

^!o G P f /<op S = N v /N c log(2/e), (46) 

and for e = 10" 5 we would have / n hl^ ^ 

\2N V /N C . This clearly shows that, if the entire self- 
energy was needed (as opposed to a few matrix elements) , 
then our proposed SGW approach would be more conve- 
nient that the standard sum-over-states approach. 



The above analysis shows that the main bottlenecks 
of our method are (i) the Fourier transform for the ap- 
plication of the Hamiltonian and (ii) the large basis sets 
adopted. In order to make the approach proposed here 
more efficient we could either move to real-space methods 
where the application of the Hamiltonian scales linearly 
with system size, 49 or reduce the size of the basis set by 
using local orbitals. 50 Fast evaluations of the operation 
Hip in order- A" operations should make it feasible GW 
calculations with A 3 t scaling and with a very favorable 
prefactor. We will come back to this aspect in Sec. VI. 

2. Green's function 

The complexity of the procedure for calculating the 
Green's function proposed in Sec. II B can be analyzed 
along the same lines of Sec. IV 1. The main differences 
in this case are that (i) the linear system Eq. (18) does 
not depend on the occupied states, (ii) the calculation is 
non self-consistent, and (iii) the calculation is performed 
for one single frequency uj + , while the entire frequency- 
dependence is generated through the multishift method. 
As a result, a factor 2NscfN v drops out of Eq. (43) and 
the computational cost of the Green's function calcula- 
tion reads: 

-^flops = 4A CG A^™. (47) 

The complexity of this calculation is significantly smaller 
than the complexity of the algorithm for the screened 
Coulomb interaction. In particular, the calculation of 
the Green's function scales as A^ t logA at . This proce- 
dure for calculating the Green's function is advantageous 
with respect to an expansion over empty states, as the 
orthogonalization of the unoccupied states would require 
a number of floating point operations scaling as ~ A 3 t . 

3. Scaling of the self-energy calculation 

The self-energy is computed in real-space after ob- 
taining G(r,r';w) and W(r, r';w) from G(G,G';a>) and 
W(G, G'; oj), respectively, and then is tranformed back 
into reciprocal space. The 6-dimensional FFT trans- 
forms require (TVq + E ) Np§ T operations for each fre- 
quency of the screened Coulomb interaction, having de- 
fined App T = 4A r SE logA r SE . The computational cost of 
this procedure scales as A at logA at , and is small with re- 
spect to the cost of calculating the screened Coulomb 
interaction. 



V. RESULTS 

In order to demonstrate the approach proposed in 
Sees. II,III we have realized a prototype implementation 
within the empirical pseudopotential method (EPM) of 
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TABLE I: Long-wavelength limit of the static sym- 
metrized inverse dielectric matrix of silicon e~^,,(q;tj) [q = 
(0.01, 0, 0)2n/a and w = 0]. We compare our calculations per- 
formed within the self-consistent Sternheimer approach with 
the results obtained in Ref. 44 using the expansion over empty 
states and the inversion of the dielectric matrix. For the cal- 
culations we sampled the Brillouin zone with 29 irreducible 
k-points, corresponding to a 8 x 8 x 8 grid, 44 ' 45 and a plane 
wave cutoff of 5 Ry. 44 Following Ref. 44 we employed the 
empirical pseudopotential parameters from Ref. 30. The re- 
ciprocal lattice vectors are in units of 2n/a, a being the lattice 
parameter. 



e Q ^,(q;w) 
G G' Ref. 44 Present work 



(0,0,0) (0,0,0) 


0.083 


0.0866 


(1,1,1) (1,1,1) 


0.605 


0.6055 


(1,1,1) (1,1,1) 


0.008 


0.0076 


(1,1,1) (1,1,1) 


0.010 


0.0102 


(1,1,1) (1,1,1) 


0.045 


0.0463 


(2,0,0) (1,1,1) 


-0.038 


-0.0382 


(2,0,0) (1,1,1) 


-0.005 


-0.0049 


(2,0,0) (2,0,0) 


0.667 


0.6671 


(2,0,0) (2,0,0) 


0.006 


0.0063 


(0,2,0) (2,0,0) 


0.016 


0.0166 



Ref. 30, and we have validated our implementation for 
the test case of silicon. 



A. Dielectric matrix 

Table I contains some of the components of the in- 
verse dielectric matrix calculated using the self-consistent 
Sternheimer method described in Sec. III. In all our 
calculations we used inverse dielectric matrices of size 
59x59, corresponding to a kinetic energy cutoff of 5 Ry 
for the screened Coulomb interaction. For the purpose of 
comparison with Ref. 44 we calculated the static and long 
wavelength limit (lu = 0, q — * 0) of the inverse dielectric 
matrix for the first few reciprocal lattice vectors. The 
authors of Ref. 44 adopted the standard approach based 
on the expansion over the unoccupied electronic states of 
the dielectric matrix, and obtained the inverse dielectric 
matrix by performing matrix inversions. In our calcula- 
tions we used the self-consistent method of Sec. Ill and no 
matrix inversion was necessary. The excellent agreement 
which can be seen in Table I between our calculations 
and Ref. 44 supports the validity of our approach. 

Next we consider the wavevector dependence of the 
head of the dielectric matrix eoo(q, (J = 0). We performed 
the calculation by using the non self-consistent method 
described in Sec. II A 2 in order to compare our results 
with Ref. 51. Figure 1 shows that our calculations are 
in very good agreement with the results of Ref. 51. The 
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0.0 0.5 1.0 1.5 2.0 

q (27i/a) [t] 

FIG. 1: (Color online) Dielectric function of silicon calculated 
using the empirical pseudopotential method and the non self- 
consistent Sternheimer method of Sec. II A 2: calculated head 
of the static dielectric function as a function of wavevector 
eoo (q, ^ = 0) (blue solid line), and results from Ref. 51 (red 
dashed line). We used a planewave kinetic energy cutoff of 5 
Ry and sampled the Brillouin zone through a uniform 8x8x8 
grid. The wavevectors are in units of 2n/a, a being the lattice 
parameter. 

slight differences between our results and those of Ref. 51 
at large wavevectors can likely be ascribed to the use 
of a limited number of empty states in the perturbative 
expansion over unoccupied states in the latter work. 

Figure 2 compares our results for the frequency depen- 
dence of the dielectric matrix at long wavelength with the 
results reported in Ref. 2. We focused in particular on 
the cases illustrated in Fig. 3 of Ref. 2. Apart from some 
small differences possibly arising from the use of the ex- 
pansion over unoccupied states in Ref. 2, even in this 
case the agreement between our calculations and those 
of Ref. 2 is very good throughout the entire frequency 
range. The agreement is consistently good for the head 
of the dielectric matrix and for diagonal and off-diagonal 
matrix elements. 



B. Pade approximants and convergence 

Figure 3 shows the quality of the analytic continuation 
from the imaginary to the real frequency axis using Pade 
approximants (cf. Appendix C). 52,53 We found that this 
procedure based on Pade approximants is generally very 
stable and requires minimal manual intervention. Ap- 
proximants of order 5 and higher are able to reproduce 
the location, the strength, and the width of the main 
plasmon-pole structure of the dielectric matrix. Head, 
wings, and body of the dielectric matrix are all described 
consistently (cf. Fig. 3). Although the singularity corre- 
sponding to the absorption onset in Fig. 3(a) is smoothed 
out by Pade approximants of low order, this effect is 
washed out when calculating the frequency convolution 
of the Green's function with the screened Coulomb inter- 
action for the GW self-energy. The advantages of per- 
forming calculations along the imaginary axis are that (i) 
the linear system in Eq. (4) becomes increasingly more 
well-conditioned when approaching large imaginary fre- 
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FIG. 2: (Color online) Frequency-dependent dielectric matrix 
of silicon at long wavelength [q = (0.01, 0, 0) 2n/a]. The cal- 
culations were performed using the empirical pseudopotential 
method of Ref. 30 and the non self-consistent Sternheimer 
method of Sec. II A 2. We used a planewave kinetic energy 
cutoff of 5 Ry and sampled the Brillouin zone thorough a 
uniform 24 x 24 x 24 grid. Such a dense Brillouin zone sam- 
pling was necessary to correctly describe the absorption on- 
set. In order to avoid null eigenvalues in the linear system 
Eq. (11) we performed the calculations by including a small 
imaginary component of 0.1 eV in the frequency U). The 
panels (a)-(c) correspond to the cases illustrated in Fig. 3 
of Ref. 2 and show eGG'(ci — ► 0, w) for (a) G = G' = 0, 
(b) G = G' = (1,1, l)27r/o, and (c) G = (1, 1, l)2n/a, 
G' = (2, 0, 0)27r/a. The blue solid lines are our calculations, 
the red dashed lines are from Ref. 2. The real and imaginary 
parts of the dielectric matrix are indicated in each panel. We 
note that the scales on the vertical axes correspond to three 
different orders of magnitude. 



quencies, and (ii) a moderate Brillouin-zone sampling is 
required to perform calculations along the imaginary axis 
unlike the case of real-axis calculations. As a result, the 
worst case scenario for the solution of the linear system 
Eq. (4) corresponds to the static case u> = 0. These tech- 
nical aspects are described in detail in Appendix B. 

The typical number of non self-consistent iterations re- 
quired to solve Eq. (4) for a fixed AVr r , w i with a relative 
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FIG. 3: (Color online) Real part of the long-wavelength di- 
electric matrix of silicon as a function of frequency. Panels (a) 
and (b) of this figure correspond to panels (a) and (b) of Fig. 2, 
respectively. The technical details of the calculations are the 
same as those described in the caption of Fig. 2. Solid blue 
lines: dielectric matrices calculated directly along the real fre- 
quency axis, from Fig. 2. Dotted, dashed, and dash-dotted 
red lines: dielectric matrices obtained from the analytic con- 
tinuation on the real axis using Pade approximants of order 6, 
11, and 26, respectively. The Pade approximants were gener- 
ated using dielectric matrices calculated along the imaginary 
frequency axis on uniform frequency grids in the range 0-50 
eV. For instance, the 6-points approximant corresponds to 
calculations at the imaginary frequencies of 0, 10, • ■ • ,50 eV. 



accuracy of £nscf = 10~ 10 using the cBiCG algorithm 
described in Appendix A is Nqg — 21 (using the precon- 
ditioner of Ref. 54). This estimate has been obtained by 
averaging over all the G, G' reciprocal lattice vectors, q- 
vectors, and imaginary frequencies. The typical number 
of self-consistent cycles required to obtain the screened 
Coulomb interaction through Eq. (4) with a relative accu- 
racy of £scf = 10~ 5 is N$cf — 5. Charge-sloshing effects 
are attenuated by using the potential mixing method pro- 
posed in Ref. 55, appropriately modified to deal with 
complex potentials. 

For completeness we report here the corresponding fig- 
ures for the calculation of the Green's function using 
Eq. (18). The average number of non self-consistent iter- 
ations required to obtain the analytic part of the Green's 
function is Nqf — 25 when preconditioning is adopted 
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(for this purpose we used a straightforward adaptation of 
the method of Ref. 54). However, multishift minimiza- 
tions as described in Appendix D do not allow for the 
use of preconditioning, and in the latter case the num- 
ber of iterations required to achieve convergence (within 
a relative accuracy £nscf = 10~ 10 ) can be as high as 
N GF ~ 100. 



C. Self-energy 

Figures 4 and 5 show the real part Re(nk|E|nk) 
and the imaginary part Im(nk|S|nk) of the GW self- 
energy calculated for the first few silicon eigenstates 
at r using our SGW method within the EPM scheme. 
Our results are compared with the calculations of of 
Ref. 18 performed within DFT/LDA and the projector- 
augmented wave method (PAW). We calculated the 
screened Coulomb interaction by using a uniform 6x6x6 
grid to sample the Brillouin zone, and Pade approximants 
of order 7 along the imaginary frequency axis. The fre- 
quency integration of the Coulomb term S c of the GW 
self-energy in Eq. (21) was performed by using a Coulomb 
cutoff ivc = 100 eV and a grid spacing of 0.5 eV. The 
Green's function was calculated using an imaginary com- 
ponent rj = 0.3 eV in Eq. (18). Apart from some differ- 
ences in the damping of the plasmaron peaks, the agree- 
ment between our calculated self-energy and the results 
of Ref. 18 is rather good throughout the entire frequency 
range ±100 eV. This finding is quite surprising since we 
are comparing our empirical pseudopotential calculations 
with low kinetic energy cutoff (5 Ry) with the ab-initio 
LDA calculations including PAW core-reconstruction of 
Ref. 18. Such agreement probably reflects the ability of 
the EPM method to provide not only a good description 
of the band structure of silicon, but also a reasonable 
description of the electronic wavefunctions. 



D. Quasi- particle excitations and spectral function 

Within the GW method the values of the quasi-particlc 
excitation energies are typically calculated by using first- 
order perturbation theory on the DFT eigenvalues. 2 The 
perturbation operator AE(r, r'; u) corresponds to the dif- 
ference between the GW self-energy E(r, r';u>) and the 
DFT exchange and correlation potential l^ xc (r). This 
approach is sensible because the complete quasi-particle 
equations are similar to the ordinary Kohn-Sham equa- 
tions if we replace the self-energy E by the DFT XC 
potential y xc . 20 

Within the EPM scheme the total potential V EPM 
acting on the electrons is specified, 30 but the electronic 
charge density is not connected to this potential through 
a self-consistent procedure. 56 This limitation makes it 
difficult to identify an XC contribution within the em- 
pirical pseudopotential. However, since the charge den- 
sity obtained within the EPM method can be regarded 
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FIG. 4: (Color online) Real part Re(nk|E(r, r'; uj)\rik) of the 
expectation value of the GW self-energy for the first 8 silicon 
eigenstates at k = 0. The solid blue lines are our SGW results 
using the empirical pseudopotential method. The dashed red 
lines are from the first-principles calculation of Ref. 18 using 
the LDA and the PAW method. Panels (a)-(d) correspond to 
the states V' 2c (Band 1), F 15c (Bands 2-4), r' 25 „ (Bands 5-7), 
and Ti„ (Band 8), respectively. The energy axis is aligned 
with the top of the valence band. The horizontal dotted lines 
indicate the calculated bare exchange contribution to the self- 
energy Re(nk|E GX (r, r'; w)|nk). 



as an approximation to the actual charge density, 57 it 
appears sensible to obtain the effective XC potential as 
a functional of the EPM charge density using the local 
density approximation. 58 ' 59 This procedure is formally 
equivalent to assuming that the unscreened ionic pseu- 
dopotential V ion is given by V" EPM - V Ha - V xc , where 
the Hartree potential V Ra and the XC potential V xc 
are calculated using the EPM charge density. This un- 
certainty on the XC potential renders the calculation of 
the quasi-particle excitation energies somewhat arbitrary, 
therefore the results presented in the following should be 
regarded as qualitative and are presented only for the pur- 
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TABLE II: GW quasi-particle excitation energies for the Y' 2c , Y 15 c, Y' 25v , Ti v states of silicon, as well as for the conduction 
band bottom of silicon close to the X\ a state. We obtained the quasi-particle energies _E„k by solving the nonlinear equations 
E n k = £ n k + R.e[AE„k(-E„k)], with A£ n k = (nk|£ — V^ c |wk). For completeness we also show the unperturbed eigenvalues 
calculated here within the EPM scheme, and those of Ref. 18 within the DFT/LDA scheme for comparison. It is interesting 
to observe that the value of the minimum band gap is not increased from the EPM value of 0.82 eV when applying the GW 
correction. In absolute terms the X\ c state shifts upwards by 0.41 eV upon the GW correction, but this shift is compensated 
by the concurrent upward shift of the Y' 25v state by 0.44 eV. This problem relates to the uncertainty on the XC potential within 
the EPM formalism, as discussed in Sec. VD. 



EPM/LDA eigenvalues Quasiparticle energies 





Present 


Ref. 18 


Present 


Ref. 2 


Ref. 18 


Expt. 


r 2c 


3.89 


3.23 


4.21 


4.08 


4.05 


4.23 a ,4.1 b 


Tl5c 


3.42 


2.54 


3.53 


3.35 


3.09 


3.40 a ,3.05 b 


^25v 


0.00 


0.00 


0.00 


0.00 


0.00 





Tlii 


-12.62 


-11.97 


-13.23 


-12.04 


-11.85 


-12.5±0.6 a 


Ale 


0.82 


0.55 


0.79 


1.29 


0.92 


1.17 a 



Ref. 60. b Ref. 61. 



pose of demonstrating a complete GW calculation within 
our SGW methodology. 

Despite the above limitations, the expectation values 
of the XC potential and of the exchange term of the 
GW self-energy calculated here are surprisingly close to 
those obtained in Ref. 2 using ab-initio pseudopotentials 
at the DFT/LDA level. Indeed, for the valence band 
top T' 2 § v state and the conduction band bottom close to 
the Xi c state our calculated XC expectation values are 
-11.27 eV and -8.97 eV, respectively while Ref. 2 gives 
-11.80 eV and -9.61 eV for the corresponding states at 
the DFT/LDA level. The agreement is even better when 
comparing the expectation values of the bare exchange 
part of the GW self-energy. In this case we find -12.43 
eV and -5.07 eV for the T' 25v state and the X lc state, 
respectively, to be compared to the corresponding values 
of -12.54 eV and -5.28 eV of Ref. 2. These results provide 
an a posteriori justification to our choice of calculating 
the XC potential using the EPM charge density and the 
LDA functional. 

Table II compares our calculated quasi-particle exci- 
tation energies for electronic states at the T point and 
for the conduction band edge of silicon with the results 
of Refs. 2,18. We find a good overall agreement between 
these different sets of calculations. Taking into account 
that we are comparing our SGW scheme within the EPM 
implementation with more sophisticated ab-initio calcu- 
lations, such agreement is rather encouraging. 

Figure 6 shows the calculated quasi-particle spectral 
function 



(nk\A\nk) = 



|ImAE 



nk | 



\uo - e„ k - ReA£ Ilk | 2 + |ImA£ nk | 2 ' 



(48) 



with A(y,v';uj) denoting the GW spectral function and 
AE„ k = (nk|E - V^ c |nk>. We note that Eq. (48) is ob- 
tained by assuming that the off-diagonal matrix elements 
of the self-energy in the basis of the unperturbed wave- 
functions are negligible. A diagonal approximation is not 



always justified, nonetheless we decided to adopt Eq. (48) 
to be consistent with Ref. 18. Figure 6 shows good overall 
agreement between our calculations and the LDA/PAW 
results of Ref. 18. This comparison demonstrates once 
again the validity of our methodology. 



VI. CONCLUSIONS AND OUTLOOK 

The results presented in Sec. V demonstrate the feasi- 
bility of the self-consistent Stcrnhcimer approach to GW 
calculations in the simple case of a prototype implemen- 
tation based on the empirical pseudopotential method. 
The extension of the present methodology to ab-initio ap- 
proaches based on norm-conserving 62 or ultrasoft 63 pseu- 
dopotentials should not present any difficulties as the cru- 
cial issues in the calculation have already been addressed 
in this work. The main advantage of the present method- 
ology consists of the definitive elimination of the unoccu- 
pied electronic states from the calculations of both the 
screened Coulomb interaction and the non-interacting 
Green's function. Another appealing aspect is that our 
methodology constitutes a generalization to frequency- 
dependent perturbations of density-functional perturba- 
tion theory, 28 which is a well established technique with 
a long history of successes. 

As discussed in Sec. IV the present approach is compa- 
rable in performance to standard 2 GW techniques. The 
question remains on whether it is possible to make sig- 
nificant improvements over the methodology proposed 
here without compromising on the numerical accuracy. 
The most time-consuming step of the entire procedure 
is the application of the single-particle Hamiltonian H 
to a search direction ip in the Hilbert space spanned by 
the wavefunction basis set during the iterative solution of 
the linear systems in Eq. (4). In order to accelerate this 
operation there are possibly three ways ahead: (i) the im- 
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FIG. 5: (Color online) Imaginary part Im(nk|E(r, r'; tj)|nk) 
of the expectation value of the GW self-energy for the first 
8 silicon eigenstates at k = 0. The solid blue lines are our 
SGW results using the empirical pseudopotential method. 
The dashed red lines are from the first-principles calculation 
of Ref. 18. Panels (a)-(d) correspond to the states Y' 2c , rise, 
P^u, and respectively. The energy axis is aligned with 
the top of the valence band. For comparison, the vertical grey 
lines indicate the locations of the logarithmic singularities at 
e„k±w p (uj p being the plasmon energy) that would arise using 
a model plasmon-pole dielectric function. 20 



provement of the minimization algorithms adopted, (ii) 
the use of sparse representations of the Hamiltonian, and 
(iii) the use of smaller basis sets. 

(i) The iterative solution of Eq. (4) is here performed 
by first solving for Aip^ at fixed AV, and then by updat- 
ing AV in the self-consistency cycle. It should be possi- 
ble, at least in principle, to combine these two operations 
in a single minimization step. This could be achieved for 
instance by using the variational formulation of density- 
functional perturbation theory developed in Ref. 64, or 
by using simulated annealing techniques such as the Car- 
Parrinello method. 65 Both of these approaches were de- 
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FIG. 6: (Color online) Expectation values of the quasiparti- 
cle spectral function (nk[ J 4(r, r'; w)|nk) for the first 8 silicon 
eigenstates at k = 0. The solid blue lines are our SGW re- 
sults using the EPM, the dashed red lines are from the first- 
principles calculation of Ref. 18 using the LDA and the PAW 
method. The expectation values have been calculated within 
the diagonal approximation of Eq. (48). In each panel the 
sharp peak near the band extrema corresponds to a well- 
defined quasi-particle, while the two broad peaks corresponds 
to plasmarons, i.e. electrons or holes coupled to a cloud of real 
plasmons. 20 The suppression of one of the plasmaron peaks 
reflects the large imaginary parts of the self-energy in the 
corresponding panels of Fig. 5. We point out the different 
vertical scale in panel (d). 



veloped for Hermitian systems, hence appropriate gener- 
alizations to non-Hcrmitian systems would be required to 
solve Eq. (4). As mentioned in Sec. VB the total number 
of cBiCG iterations required to solve Eq. (4) is typically 
-^CG-^SCF ~ 100, therefore the variational formulation 
of DFPT or the Car-Parrinello minimization would be 
convenient if they resulted in a significant reduction of 
the number of iterations over this figure. 

(ii) Another possibility for improving the procedure 
presented in this work is to resort to a real-space rep- 
resentation of the kinetic energy operator in the single- 
particle Hamiltonian. This can be achieved by calculat- 
ing the kinetic energy on a real-space mesh using finite- 
differences methods. 49 The main advantage of this ap- 
proach would be to have GW calculations scaling with 
the cube of the system size N^ t instead of N^ t \ogN at - 
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However, in actual calculations the numerical prefactor 
associated with this scaling could be unfavorable. Indeed, 
the fast Fourier transforms used in a planewaves repre- 
sentation to calculate the product of the potential and 
the wavefunctions require 2 • 4A^ icn logA^ icn floating point 
operations, while the cost of the real-space calculation of 
the Laplacian operator to sixth order is 37iV^ cn (for a 
supercell with orthogonal axes). 49 While these estimates 
seem to speak in favor of the real-space method, it is 
advisable to consider that the preconditioner used here 54 
for the cBiCG minimization cannot be simply adapted to 
real-space calculations. In absence of effective precondi- 
tioners the planewave method remains advantageous for 
essentially any relevant system size. 

(iii) Another interesting option for improving our 
methodology is to drastically reduce the size of the ba- 
sis set adopted. This could be achieved for instance 
by adapting our implementation to electronic structure 
packages exploiting local orbitals basis sets. 50 Interest- 
ingly the three possibilities here outlined are not mutu- 
ally exclusive, and probably an appropriate combination 
of all of these would eventually open the way to the study 
of electronic excitations in very large systems using the 
GW method. 

In summary, we propose a new methodology for per- 
forming GW calculations using the self-consistent Stern- 
heimer equation (SGW). We show how to calculate the 
screened Coulomb interaction and the non-interacting 
Green's function without resorting to unoccupied elec- 
tronic states. We successfully demonstrate our method 
within a planewaves empirical pseudopotential imple- 
mentation and compare with previous studies for the 
prototypical test case of silicon. In our method the stan- 
dard generalzied plasmon-pole approximation for the fre- 
quency dependence of the screened Coulomb interaction 
has been replaced by a direct calculation along the imag- 
inary frequency axis, followed by an analytic continua- 
tion to the real axis. In addition, we introduce the use 
of multishift linear system solvers for the simultaneous 
calculation of multiple frequency responses at the cost of 
one single iterative minimization. 

It is our plan to adapt the present approach to deal 
with first-principles pseudopotentials, and to explore the 
performance of our procedure in a local orbital real-space 
implementation. We hope that this work will stimulate 
further effort to develop improved methodologies for ex- 
cited states calculations in large systems. 
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APPENDIX A: PRECONDITIONED COMPLEX 
BICONJUGATE GRADIENT METHOD 

In this work we solve the linear systems in Eq. (4) 
using the complex biorthogonal variant (cBiCG) of the 
conjugate gradient method (CG) following Ref. 47. This 
method is an extension of the standard conjugate gradi- 
ents algorithm to the case of general complex matrices. 
As the cBiCG algorithm introduced in Ref. 47 did not in- 
clude preconditioning, in this appendix we describe the 
preconditioned version of the algorithm which we derived 
following Ref. 66. We are interested in solving the linear 
system 

Ax = b, (Al) 

with A a complex linear operator (not necessarily Her- 
mitian), b a complex vector, and x the unknown solution 
vector. In the cBiCG algorithm of Ref. 47 two sequences 
of residuals r n and f„ are generated in such a way that 
successive residuals are biorthogonal [i.e. (r„ + i|f„) = 
and {f n +i\r n ) = 0]. Two sequences of search directions 
p n and p n are generated so that successive directions are 
biconjugate [i.e. (Ap n+1 \p n ) = and (A^p n+1 \p n ) = 0]. 
The algorithm starts by setting the initial residuals to 
ro = b — Axq (xq being the initial guess for the solution 
vector x) and f — t-q, and the initial search directions 
to po = r and p = Pq. Next the solution vector, the 
search directions, and the residuals are updated at each 



iteration as follows: 

a n = (f n \r n ) I (p n \Ap n ) (A2) 

x n +i = x n + a n p n (A3) 

r n +i = r n - a n Ap n (A4) 

f n +i = f n - a* n A^p n (A5) 

p n = -(A^p n \r n+1 )/(p n \A Pn ) (A6) 

Pn+l = r n+ i + (3 n p n (Al) 

p n+ i = f n+1 +(3*p n . (A8) 



The time-consuming step in this algorithm corresponds 
to the application of the operators A and to the search 
directions p n and p n . As there are two such operations 
per iteration, the computational complexity is twice that 
of the standard conjugate gradient algorithm. 

The preconditioning of the linear operator can be 
achieved by left-multiplying the linear system in Eq. (Al) 
by M^ 1 : M~ x Ax = M~ x b. If we assume that the pre- 
conditioner M can be written as M — E T E, then we can 
rewrite the system as follows: 

E- 1 AE- T E T x = E~H. (A9) 
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By defining A' = E' 1 AE~ T , x' = E T x, and b' = E^b 
we obtain the transformed system A'x' = b', for which 
the standard cBiCG method applies. While this proce- 
dure is formally correct, it is not convenient to explicitly 
transform the linear operator. It is convenient instead 
to rewrite the procedure in terms of A, b, and x by per- 
forming a few formal manipulations. For this purpose we 
make the substitutions r' = E~ 1 r and p' = E T p. Some 
algebra leads straightforwardly to the preconditioned ver- 
sion of the cBiCG algorithm: 



On = (rnlM-^/iPnlApn) (A10) 

x n+1 = x n +a n p n (All) 

r„+i = r n - a n Ap n (A12) 

f n +i = f n - a* n A^p n (A13) 

(3 n = -(AtpnlM-Vn+i)/^!^) (A14) 

Pn+1 = M~ 1 r n+1 + /3 n p n (A15) 

p n+1 = M-^n+i + P*p n . (A16) 

The preconditioned cBiCG algorithm needs to be intial- 

ized with ro = b — Axq, po = M ro, fo = r^, and 



Po = Po- In this work we have used the Teter-Payne- 
Allan function as the preconditioner M -1 . 54 



APPENDIX B: CONDITION NUMBER OF THE 
LINEAR SYSTEM 

a. Screened Coulomb interaction 

The iterative calculation of the screened Coulomb in- 
teraction through Eq. (4) at finite real frequencies lj can 
be considerably more time-consuming than in the static 
(cj = 0) case. Simple tests indicate that the number of 
iterations required to achieve convergence increases with 
increasing frequency lj. This behavior suggests that the 
linear operator in Eq. (4) becomes progressively more ill- 
conditioned as the frequency to increases. 

In order to rationalize this observation, we here exam- 
ine the condition number of the linear operator H—e v ±u> 
in Eq. (4). The minimum number of iterations N m { n re- 
quired for the solution of a linear system using the con- 
jugate gradients algorithm is given by 

iVmin - ^log(2/e), (Bl) 

k being the condition number of the linear operator and 
£ the desired relative accuracy. 66 While the estimate 
Eq. (Bl) has been derived for the original CG algorithm, 
we found empirically that it also provides a reasonable 
description of the convergence rate of the complex cBiCG 
version. The condition number k of a linear operator can 
be calculated as the ratio of its largest to smallest eigen- 
values. For a given valence state \v') the linear opera- 
tor H — e v > + aP occ — u> in Eq. (7) has the eigenvalues 
s v — e v > + a — uj and e c — e v > — u>. 



Let consider first the simplest case where u> = and 
a > 0. In this case we find by inspection that the small- 
est eigenvalue is mm(E g , \ a — W OC c|)> E g being the fun- 
damental energy gap and W occ the valence bandwidth. 
It is common practice to set a = 2W occ to avoid null 
eigenvalues. 28 With this choice the smallest eigenvalue 
becomes E g . On the other hand, the largest eigenvalue 
can be approximated by the cutoff energy of the wave- 
function basis set E cut . In this case the condition num- 
ber reads k = E cut /E g . As an example, if we are using 
a plane-waves basis with a kinetic enegry cutoff of 30 
Ry, we have an electron energy gap of 1 cV, and the 
our desired accuracy is e = 10~ 10 , then according to 
Eq. (Bl) the minimum number of iterations required to 
solve the linear system would be iV m ; n = 240. Empirical 
tests show that this estimate is quite accurate for the sys- 
tem considered in the present work. In order to improve 
the convergence rate it is useful to employ precondition- 
ing techniques. We here adopt the Teter-Payne- Allan 
preconditioner 54 in order to "compress" the eigenvalue 
spectrum and thereby reduce the condition number. An 
ideal preconditioner would make the linear operator per- 
fectly well conditioned (k = 1). In this case the optimal 
number of iterations (for a relative accuracy e = 10 -10 ) 
would be as small as -/V m i n . pc = 12. We have found empir- 
ically that by using the Teter-Payne- Allan preconditioner 
the number of iterations required to achieve convergence 
was in all cases in the range A^tpa =15-40. 

We now consider the case of u> < 0. Simple algebra 
shows that in this case k = E cxA /(E g + w) when a = 
2 Wo CC . Hence in this case the larger the frequency u>, the 
better conditioned the linear system. We checked this 
result by explicit calculations. 

The worst case in terms of condition number is found 
when uj > 0. In fact, as soon as the frequency exceeds the 
optical excitation threshold u > E g , the linear operator 
acquires null eigenvalues corresponding to the resonance 
condition lj = e c — e v > . In this latter case the condition 
number k(u>) exhibits significant structure, reflecting the 
joint density of states of the system. Even after precondi- 
tioning the system, the number of iterations required to 
achieve convergence can be as high as A m ; n = 500, thus 
rendering this avenue unpractical. The calculation of the 
screened Coulomb interaction for frequencies slightly off 
the real axis u) + iff, with lj > and small 77, leads to 
only a small improvement of the convergence rate. The 
difficulty of solving iteratively the linear system Eq. (4) 
for large positive frequencies is accompanied by the ad- 
ditional difficulty of adequately sampling the Brillouin 
zone to describe the singularities at lj = e c — e v > . 

Altogether these considerations suggest that an itera- 
tive solution of the linear system along the real axis is 
not convenient from the computational viewpoint. For 
this reason we decided to evaluate the screened Coulomb 
interaction along the imaginary axis and then to analyt- 
ically continue the functions to the real axis using Pade 
approximants. 18 ' 52,53 The motivation behind our choice 
becomes obvious after considering a simple plasmon-pole 
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model of the screened Coulomb interaction: 2 



v + 



W -v 



L0 V 



(B2) 



* lid + LO p lO — lO p . 

where io p is the pole frequency and Wo the static screened 
Coulomb interation. Analytical continuation of this func- 
tion to the imaginary axis yields 

W -v 



W(u =i/3) = v + 



1 + ((3/u> p ) 



2 ' 



(B3) 



Equation (B3) indicates that the screened Coulomb in- 
teraction along the imaginary axis contains the same 
amount of information as the one on the real axis (no p 
and Wo), and at the same time does not exhibit any 
singularities. In this case the condition number reads 
(assuming no preconditioning and a = for simplicity) 
k = [(El + (3 2 )/(E 2 ut + /3 2 )]3, and tends to unity for 
large imaginary frequencies. As a result, the worst case 
scenario for the solution of the linear system corresponds 
to the static case u = 0. 

In summary by solving iteratively the linear system 
along the imaginary axis we circumvent the difficulties 
associated with the ill-conditioning of the linear system in 
Eq. (4) occurring at real frequencies and the necessity of 
dense Brillouin zone sampling. The details of the analytic 
continuation are discussed in Appendix C. 



b. Green's function 

A similar analysis can be carried out for the calcu- 
lation of the Green's function using the method intro- 
duced in Sec. II B. It is straightforward to establish that 
in this case the condition number of the system is given 
by k = E cut /5. As the infinitesimal 6 is typically taken 
to be 0.1 eV, we are effectively dealing with a situation 
analogous to a small-bandgap semiconductor. The TPA 
preconditioner can be adopted to reduce the condition 
number to k = E^P M /6, where -E^„ M is the expectation 
value of the kinetic energy of the highest occupied state 
and is independent of the basis set cutoff. Numerical 
tests confirm that this is indeed a sensible and effective 
strategy. 



APPENDIX C: ANALYTIC CONTINUATION 
USING PADE APPROXIMANTS 

In order to perform the analytic continuation of 
the screened Coulomb interaction from the imagi- 
nary axis to the real axis, we employ diagonal Pade 
approximants. 18 ' 52,53 The Pade approximant of order N 
is the optimal rational approximation to a target func- 
tion f(oj) known in N distinct points w n , n = 1, • • • , N. 
When N is an odd integer the diagonal Pade approxi- 
mant reads 

, , Po+piu + ---+p {N _ 1)/2 J N - 1 )/ 2 

= - ; : ; — . .(at-i)/2 ' ( C1 ) 



1 + quo H h 9(jv-i)/2W 



and its coefficients p ,Pi,--- ,f>(;v-i)/2, <7i,-" ,9(jv-i)/2 
are determined by matching the approximant to the tar- 
get function in N points P/v(w„) = f(to n ), n = 1, • • • ,N. 
Both the coefficients and the Pade approximant can 
be calculated very efficiently using a simple recursive 
algorithm. 53 Some experimentation indicates that ap- 
proximants of order N > 5 are necessary to reproduce a 
plasmon-polc spectral shape including a finite lincwidth. 
This observation can be rationalized by considering that 
a plasmon-pole spectral function is completely defined 
by the values of the function at w = 0, the location, 
strength, and width of the pole, and the asymptotic value 
at to = +ioo. Some of this information is redundant and 
can be obtained by using sum rules. 2 The parity of the 
screened Coulomb interaction can also be exploited to 
minimize the number of input frequencies. The advan- 
tage of the Pade approximant is that a more refined de- 
scription of the frequency-dependent screened Coulomb 
interaction can simply be achieved by calculating addi- 
tional points along the imaginary axis. 

We also investigated the possibility of analytically con- 
tinuing the screened Coulomb interaction by using a 
multi-pole expansion as suggested in Ref. 17. We tried 
one-, two-, and three-pole expansions by determining 
the coefficients using the simplex method of Nclder and 
Mead. 67 The single-pole approximation appears robust 
but the quality of the real-axis continuation is poorer 
than what we obtained by using Pade approximants. 
Multi-pole approximations were found to be unreliable 
because of their high sensitivity to the initial guesses 
for the coefficients. Our experience therefore is that the 
multi-pole expansion is not an optimal choice for an au- 
tomated procedure where the analytic continuation has 
to be performed for every G, G', and q of the screened 
Coulomb interaction without manual intervention. 



APPENDIX D: SIMULTANEOUS CALCULATION 
OF THE SUSCEPTIBILITY AT MULTIPLE 
FREQUENCIES 

The linear systems Eqs. (11) and (18) can be solved 
efficiently by using the "multishift" cBiCG method of 
Ref. 29. Multishift methods exploit the knowledge gained 
during the iterative solution of the seed system Ax = b 
to determine the solutions of the shifted system Ax + 
lux = b with only a small computational overhead. The 
rationale behind such method is that the seed system 
and the shifted system share the same Krilov subspaces 
{b, Ab, A 2 b, ■ ■ ■ }, therefore the residuals of the seed and 
of the shifted systems can be taken to be collinear. 29 

This multishift technique allows us to determine the 
entire frequency-dependence of the dielectric matrix by 
performing one single static calculation for each set of 
parameters [q,G] in Eq. (34). For the seed system the 
algorithm is still given by Eqs. (A2)-(A8). For the shifted 
system we replace the calculation of the residuals r n>u 
and of the coefficients a n _ u , f3 n , u corresponding to the 
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frequency u> by the following relations: 



(Dl) 

where the scaling factor Tr n +i, u is calculated through the 
recurrence relation 

= (l+wa„)7r„. w + — ^ — (7r„, w -7r n _i,a,). (D2) 

Q!n-1 

In order to obtain collinear residuals r n and r n ,ui we need 
to initialize the algorithm using xo = 0. 

The use of Eqs. (Dl) and (D2) allows us to skip the 
time-consuming operations involving the Hamiltonian in 
Eqs. (A2) and (A6). This method is extremely conve- 
nient for determining the frequency-dependent suscepti- 
bility for many frequencies at the cost of one single cal- 
culation. 

We point out that this method still carries some draw- 
backs. One limitation is that this method cannot be ap- 



plied to the self-consistent system of Eq. (4), because the 
know-term on the right-hand side depends on the fre- 
quency uj itself. Therefore the use of the shifted cBiCG 
method is only possible for non self- consistent calcula- 
tions of the dielectric matrix and requires explicit matrix 
inversions to determine the screened Coulomb interac- 
tion. This approach can be regarded as an improved 
version of the technique proposed in Ref. 24. 



Another limitation is that the shifted cBiCG method 
does not allow for the use of preconditioners. In fact 
the preconditioned seed system M~ x Ax = M~ x b and the 
preconditioned shifted system M~ 1 Ax+u)M~ 1 x = M~ x b 
do not share the same Krilov subspaces, hence the resid- 
uals cannot be taken to be collinear. 68 The practical con- 
sequence is that for systems with large basis set energy 
cutoffs and small band gaps, the number of iterations 
required to achieve convergence could be impractically 
large. 
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